TFG JOANA
Carregar llibreries
Carregar la base de dades original (corregida) i eliminar les dades de Mordèl·lids i de les zones “Llig”, “Mura”, i “TV-7041”
# Obtenir el directori de l'arxiu actual
directorio <- getwd()
dades_finals <- read_excel(file.path(directorio,"Carreteres_ITS+altres.xlsx"))
dades_finals <- dades_finals %>%
filter(Familia != "Mordellidae" & !Zona %in% c("Llig", "Mura", "TV-7041"))
head(dades_finals)## # A tibble: 6 × 25
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 28 Campus Juliol 7 2021 4 NS <NA>
## # ℹ 17 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>
ANÀLISI DESCRIPTIVA
# Definir nova variable: "Mostreig", segons Zona, Mes i Any
dades_finals <- dades_finals %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-"))Riquesa
Riquesa total
## # A tibble: 1 × 1
## riquesa_total
## <int>
## 1 137
Riquesa per cada mostreig i tractament
## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
## # A tibble: 11 × 3
## # Groups: Mostreig [6]
## Mostreig Tractament riquesa
## <chr> <chr> <int>
## 1 Campus-Juliol-2021 NS 20
## 2 Campus-Juliol-2021 S 21
## 3 Campus-Jun-2021 NS 45
## 4 Campus-Jun-2021 S 34
## 5 Campus-juliol-2020 NS 39
## 6 Campus-juliol-2020 S 24
## 7 Campus-maig-2020 NS 45
## 8 Franqueses-Maig-2021 NS 20
## 9 Franqueses-Maig-2021 S 11
## 10 Sabadell-Juny-2021 NS 16
## 11 Sabadell-Juny-2021 S 12
Riquesa per cada zona i tractament
riquesa_zones <- dades_finals %>%
group_by(Zona, Tractament) %>%
summarise(riquesa = n_distinct(ID))## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 3
## # Groups: Zona [3]
## Zona Tractament riquesa
## <chr> <chr> <int>
## 1 Campus NS 98
## 2 Campus S 59
## 3 Franqueses NS 20
## 4 Franqueses S 11
## 5 Sabadell NS 16
## 6 Sabadell S 12
Riquesa per grups
# Abelles
riquesa_abelles<- dades_finals %>%
filter(Ordre == 'HymenopteraAbella') %>%
summarise(riquesa_abelles=n_distinct(ID))
riquesa_abelles## # A tibble: 1 × 1
## riquesa_abelles
## <int>
## 1 80
# Coleòpters
riquesa_coleoptera<- dades_finals %>%
filter(Ordre == 'Coleoptera') %>%
summarise(riquesa_coleoptera=n_distinct(ID))
riquesa_coleoptera## # A tibble: 1 × 1
## riquesa_coleoptera
## <int>
## 1 33
# Sírfids
riquesa_sirfids<- dades_finals %>%
filter(Familia == 'Syrphidae') %>%
summarise(riquesa_sirfids=n_distinct(ID))
riquesa_sirfids## # A tibble: 1 × 1
## riquesa_sirfids
## <int>
## 1 9
Abundància
Abundància total
## # A tibble: 1 × 1
## Numero_total
## <int>
## 1 1111
Abundància per cada zona i tractament
## `summarise()` has grouped output by 'Zona'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 3
## # Groups: Zona [3]
## Zona Tractament Numero_total
## <chr> <chr> <int>
## 1 Campus NS 649
## 2 Campus S 323
## 3 Franqueses NS 45
## 4 Franqueses S 14
## 5 Sabadell NS 29
## 6 Sabadell S 51
Abundància de cada gènere total
## # A tibble: 6 × 2
## Genere Numero_total
## <chr> <int>
## 1 Acmaeodera 6
## 2 Amegilla 1
## 3 Andrena 6
## 4 Anthaxia 29
## 5 Anthidium 10
## 6 Apis 11
Abundància de cada ID per zona i tractament
abundaciazona_ID <- dades_finals %>%
group_by(Zona, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Tractament'. You can override
## using the `.groups` argument.
## # A tibble: 6 × 4
## # Groups: Zona, Tractament [1]
## Zona Tractament ID Numero_total
## <chr> <chr> <chr> <int>
## 1 Campus NS Acmaeodera cylindrica 5
## 2 Campus NS Amegilla sp. 1
## 3 Campus NS Anthaxia (Anthaxia) 1
## 4 Campus NS Anthaxia (Melanthaxia) 5
## 5 Campus NS Anthaxia hungarica 10
## 6 Campus NS Anthaxia sp. 7
Abundancia de cada ID per cada mostreig i tractament
abundacia_ID <- dades_finals %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
## # A tibble: 6 × 4
## # Groups: Mostreig, Tractament [1]
## Mostreig Tractament ID Numero_total
## <chr> <chr> <chr> <int>
## 1 Campus-Juliol-2021 NS Amegilla sp. 1
## 2 Campus-Juliol-2021 NS Anthaxia sp. 1
## 3 Campus-Juliol-2021 NS Bombus pascuorum 1
## 4 Campus-Juliol-2021 NS Ceratina cucurbitina 1
## 5 Campus-Juliol-2021 NS Ceratina parvula 7
## 6 Campus-Juliol-2021 NS Chrysididae 1
DIVERSITAT DE SHANNON I CORBES D’ACUMULACIÓ D’ESPÈCIES
Diversitat de Shannon per Zona i Tractament
# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons la zona i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundaciazona_ID, Zona + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0# Agrupar Zona amb Tractament (en una nova variable Agrupacio2)
data_shannon_ZT <- taula_aux2 %>%
mutate(Agrupacio2 = paste(Zona, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Zona, -Tractament)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio2 i establir-la com a identificador de les files
rownames(data_shannon_ZT) <- data_shannon_ZT$Agrupacio2
data_shannon_ZT <- data_shannon_ZT %>% dplyr::select(-Agrupacio2) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv2 <- diversity(data_shannon_ZT)
print(shannondiv2)## Campus_NS Campus_S Franqueses_NS Franqueses_S Sabadell_NS
## 3.612943 2.870658 2.627033 2.341994 2.382088
## Sabadell_S
## 2.003690
Diversitat de Shannon per Mostreig i Tractament
# Preparar dades: Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0# Agrupar Mostreig amb Tractament (en una nova variable Agrupacio)
data_shannon_MT <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MT) <- data_shannon_MT$Agrupacio
data_shannon_MT <- data_shannon_MT %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig i Tractament
shannondiv <- diversity(data_shannon_MT)
print(shannondiv)## Campus-juliol-2020_NS Campus-juliol-2020_S Campus-Juliol-2021_NS
## 3.011325 2.404801 2.593731
## Campus-Juliol-2021_S Campus-Jun-2021_NS Campus-Jun-2021_S
## 1.933200 3.161179 2.856531
## Campus-maig-2020_NS Franqueses-Maig-2021_NS Franqueses-Maig-2021_S
## 2.957046 2.627033 2.341994
## Sabadell-Juny-2021_NS Sabadell-Juny-2021_S
## 2.382088 2.003690
Corba d’acumulació d’espècies per mostreig (11 mostrejos)
## Warning in cor(x > 0): the standard deviation is zero
Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 137 236.8148 31.40658 207 26.07681 250.1182 167.3704 12.88879 11
En aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 236 (chao), 207 (jackknife primer ordre), 250 (jackknife segon ordre), i 167 (bootstrap).
Diversitat de Shannon per Mostreig, Tractament i Trampa
# Preparar dades:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa
abundacia_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
## # A tibble: 564 × 5
## # Groups: Mostreig, Tractament, Trampa [101]
## Mostreig Tractament Trampa ID Numero_total
## <chr> <chr> <dbl> <chr> <int>
## 1 Campus-Juliol-2021 NS 2 Lasioglossum glabriusculum 1
## 2 Campus-Juliol-2021 NS 4 Curculionidae 1
## 3 Campus-Juliol-2021 NS 4 Lasioglossum malachurum 1
## 4 Campus-Juliol-2021 NS 4 Pompilidae 1
## 5 Campus-Juliol-2021 NS 4 Pyronia cecilia 1
## 6 Campus-Juliol-2021 NS 6 Lasioglossum glabriusculum 1
## 7 Campus-Juliol-2021 NS 6 Lasioglossum politum 3
## 8 Campus-Juliol-2021 NS 6 Pompilidae 1
## 9 Campus-Juliol-2021 NS 8 Chrysididae 1
## 10 Campus-Juliol-2021 NS 8 Hylaeus sp. 1
## # ℹ 554 more rows
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundacia_trampa, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0# Agrupar Mostreig amb Tractament i Trampa
data_shannon_MTTR <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_MTTR) <- data_shannon_MTTR$Agrupacio
data_shannon_MTTR <- data_shannon_MTTR %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondiv <- diversity(data_shannon_MTTR)
print(shannondiv)## Campus-juliol-2020_NS_2 Campus-juliol-2020_NS_3 Campus-juliol-2020_NS_4
## 1.7917595 1.8310205 2.0449312
## Campus-juliol-2020_NS_6 Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10
## 1.7328680 1.8310205 1.0735428
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13
## 0.0000000 1.7478681 1.6434177
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16
## 1.1240357 1.5498260 1.3862944
## Campus-juliol-2020_NS_20 Campus-juliol-2020_S_1 Campus-juliol-2020_S_5
## 2.2718685 0.6931472 1.9072840
## Campus-juliol-2020_S_7 Campus-juliol-2020_S_9 Campus-juliol-2020_S_17
## 0.6924212 0.0000000 1.0986123
## Campus-juliol-2020_S_18 Campus-juliol-2020_S_19 Campus-juliol-2020_S_21
## 1.6094379 1.5607104 1.1189689
## Campus-juliol-2020_S_22 Campus-Juliol-2021_NS_2 Campus-Juliol-2021_NS_4
## 1.0751393 0.0000000 1.3862944
## Campus-Juliol-2021_NS_6 Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10
## 0.9502705 1.4648164 0.6931472
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16
## 1.3296613 1.3321790 0.0000000
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20 Campus-Juliol-2021_S_1
## 1.3862944 1.0986123 1.0986123
## Campus-Juliol-2021_S_3 Campus-Juliol-2021_S_5 Campus-Juliol-2021_S_7
## 0.6365142 2.0963526 1.5000764
## Campus-Juliol-2021_S_9 Campus-Juliol-2021_S_11 Campus-Juliol-2021_S_13
## 1.1973401 1.4388830 0.0000000
## Campus-Juliol-2021_S_15 Campus-Juliol-2021_S_17 Campus-Jun-2021_NS_2
## 0.6931472 0.5091373 2.3345491
## Campus-Jun-2021_NS_4 Campus-Jun-2021_NS_6 Campus-Jun-2021_NS_8
## 2.0854684 1.1537419 2.6197294
## Campus-Jun-2021_NS_10 Campus-Jun-2021_NS_12 Campus-Jun-2021_NS_14
## 1.8310205 2.5575077 1.7623137
## Campus-Jun-2021_NS_16 Campus-Jun-2021_NS_18 Campus-Jun-2021_NS_20
## 1.3593685 1.9061547 0.9502705
## Campus-Jun-2021_S_1 Campus-Jun-2021_S_3 Campus-Jun-2021_S_5
## 2.2718685 2.0028830 1.7351265
## Campus-Jun-2021_S_7 Campus-Jun-2021_S_9 Campus-Jun-2021_S_11
## 0.9502705 2.1383972 1.5595812
## Campus-Jun-2021_S_13 Campus-Jun-2021_S_15 Campus-Jun-2021_S_17
## 1.8576598 0.7356219 2.1639557
## Campus-Jun-2021_S_19 Campus-maig-2020_NS_1 Campus-maig-2020_NS_2
## 1.9722470 0.8392696 1.3642812
## Campus-maig-2020_NS_3 Campus-maig-2020_NS_4 Campus-maig-2020_NS_5
## 1.6957425 2.5521719 1.5171064
## Campus-maig-2020_NS_6 Campus-maig-2020_NS_7 Campus-maig-2020_NS_8
## 2.2614951 1.2424533 1.7480673
## Campus-maig-2020_NS_9 Campus-maig-2020_NS_10 Campus-maig-2020_NS_11
## 1.6769878 1.9722470 1.7480673
## Campus-maig-2020_NS_12 Campus-maig-2020_NS_13 Campus-maig-2020_NS_14
## 1.0986123 1.5481155 1.8891592
## Campus-maig-2020_NS_15 Campus-maig-2020_NS_16 Campus-maig-2020_NS_17
## 1.8392967 1.0356097 1.3296613
## Campus-maig-2020_NS_18 Campus-maig-2020_NS_19 Campus-maig-2020_NS_20
## 1.5810938 1.3535914 0.9002561
## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3
## 1.0789922 0.8675632 0.6931472
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5 Franqueses-Maig-2021_S_1
## 2.2739657 0.7549968 1.3321790
## Franqueses-Maig-2021_S_2 Franqueses-Maig-2021_S_4 Sabadell-Juny-2021_NS_1
## 1.0986123 1.5607104 1.7478681
## Sabadell-Juny-2021_NS_4 Sabadell-Juny-2021_NS_7 Sabadell-Juny-2021_NS_8
## 0.9502705 1.7917595 1.4750763
## Sabadell-Juny-2021_NS_11 Sabadell-Juny-2021_S_2 Sabadell-Juny-2021_S_3
## 1.3862944 0.6365142 0.5623351
## Sabadell-Juny-2021_S_5 Sabadell-Juny-2021_S_6 Sabadell-Juny-2021_S_9
## 1.4995094 1.0986123 0.6931472
## Sabadell-Juny-2021_S_10 Sabadell-Juny-2021_S_12
## 0.5004024 1.2852930
Corba d’acumulació d’espècies segons mostreig, tractament i trampa
Extrapolació de dades (quantes espècies més trobaríem si féssim més mostrejos)
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 137 224.1493 29.51259 201.3564 11.27322 241.7798 164.5183 5.987207 101
Amb aquest mostreig s’han trobat 137 espècies (riquesa), però s’estima que hi ha altres espècies “unseen” (no detectades). Segons diferents estimacions, la riquesa total tenint en compte aquestes espècies no detectades seria 224 (chao), 201 (jackknife primer ordre), 241 (jackknife segon ordre), i 164 (bootstrap).
Corba i extrapolació pel Campus:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades pel Campus
abundacia_trampa_campus <- dades_finals %>%
filter(Zona == 'Campus') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxcampus <- cast(abundacia_trampa_campus, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxcampus <- as.data.frame(taula_auxcampus)
taula_auxcampus[is.na(taula_auxcampus)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_campus <- taula_auxcampus %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_campus) <- data_shannon_campus$Agrupacio
data_shannon_campus <- data_shannon_campus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivcampus <- diversity(data_shannon_campus)
print(shannondivcampus)## Campus-juliol-2020_NS_2 Campus-juliol-2020_NS_3 Campus-juliol-2020_NS_4
## 1.7917595 1.8310205 2.0449312
## Campus-juliol-2020_NS_6 Campus-juliol-2020_NS_8 Campus-juliol-2020_NS_10
## 1.7328680 1.8310205 1.0735428
## Campus-juliol-2020_NS_11 Campus-juliol-2020_NS_12 Campus-juliol-2020_NS_13
## 0.0000000 1.7478681 1.6434177
## Campus-juliol-2020_NS_14 Campus-juliol-2020_NS_15 Campus-juliol-2020_NS_16
## 1.1240357 1.5498260 1.3862944
## Campus-juliol-2020_NS_20 Campus-juliol-2020_S_1 Campus-juliol-2020_S_5
## 2.2718685 0.6931472 1.9072840
## Campus-juliol-2020_S_7 Campus-juliol-2020_S_9 Campus-juliol-2020_S_17
## 0.6924212 0.0000000 1.0986123
## Campus-juliol-2020_S_18 Campus-juliol-2020_S_19 Campus-juliol-2020_S_21
## 1.6094379 1.5607104 1.1189689
## Campus-juliol-2020_S_22 Campus-Juliol-2021_NS_2 Campus-Juliol-2021_NS_4
## 1.0751393 0.0000000 1.3862944
## Campus-Juliol-2021_NS_6 Campus-Juliol-2021_NS_8 Campus-Juliol-2021_NS_10
## 0.9502705 1.4648164 0.6931472
## Campus-Juliol-2021_NS_12 Campus-Juliol-2021_NS_14 Campus-Juliol-2021_NS_16
## 1.3296613 1.3321790 0.0000000
## Campus-Juliol-2021_NS_18 Campus-Juliol-2021_NS_20 Campus-Juliol-2021_S_1
## 1.3862944 1.0986123 1.0986123
## Campus-Juliol-2021_S_3 Campus-Juliol-2021_S_5 Campus-Juliol-2021_S_7
## 0.6365142 2.0963526 1.5000764
## Campus-Juliol-2021_S_9 Campus-Juliol-2021_S_11 Campus-Juliol-2021_S_13
## 1.1973401 1.4388830 0.0000000
## Campus-Juliol-2021_S_15 Campus-Juliol-2021_S_17 Campus-Jun-2021_NS_2
## 0.6931472 0.5091373 2.3345491
## Campus-Jun-2021_NS_4 Campus-Jun-2021_NS_6 Campus-Jun-2021_NS_8
## 2.0854684 1.1537419 2.6197294
## Campus-Jun-2021_NS_10 Campus-Jun-2021_NS_12 Campus-Jun-2021_NS_14
## 1.8310205 2.5575077 1.7623137
## Campus-Jun-2021_NS_16 Campus-Jun-2021_NS_18 Campus-Jun-2021_NS_20
## 1.3593685 1.9061547 0.9502705
## Campus-Jun-2021_S_1 Campus-Jun-2021_S_3 Campus-Jun-2021_S_5
## 2.2718685 2.0028830 1.7351265
## Campus-Jun-2021_S_7 Campus-Jun-2021_S_9 Campus-Jun-2021_S_11
## 0.9502705 2.1383972 1.5595812
## Campus-Jun-2021_S_13 Campus-Jun-2021_S_15 Campus-Jun-2021_S_17
## 1.8576598 0.7356219 2.1639557
## Campus-Jun-2021_S_19 Campus-maig-2020_NS_1 Campus-maig-2020_NS_2
## 1.9722470 0.8392696 1.3642812
## Campus-maig-2020_NS_3 Campus-maig-2020_NS_4 Campus-maig-2020_NS_5
## 1.6957425 2.5521719 1.5171064
## Campus-maig-2020_NS_6 Campus-maig-2020_NS_7 Campus-maig-2020_NS_8
## 2.2614951 1.2424533 1.7480673
## Campus-maig-2020_NS_9 Campus-maig-2020_NS_10 Campus-maig-2020_NS_11
## 1.6769878 1.9722470 1.7480673
## Campus-maig-2020_NS_12 Campus-maig-2020_NS_13 Campus-maig-2020_NS_14
## 1.0986123 1.5481155 1.8891592
## Campus-maig-2020_NS_15 Campus-maig-2020_NS_16 Campus-maig-2020_NS_17
## 1.8392967 1.0356097 1.3296613
## Campus-maig-2020_NS_18 Campus-maig-2020_NS_19 Campus-maig-2020_NS_20
## 1.5810938 1.3535914 0.9002561
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 118 184.0553 24.55884 170.3457 9.341038 201.8116 140.7347 5.098803 81
Corba i extrapolació per Franqueses:
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Franqueses
abundacia_trampa_franqueses <- dades_finals %>%
filter(Zona == 'Franqueses') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxfranq <- cast(abundacia_trampa_franqueses, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxfranq <- as.data.frame(taula_auxfranq)
taula_auxfranq[is.na(taula_auxfranq)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_franqueses <- taula_auxfranq %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_franqueses) <- data_shannon_franqueses$Agrupacio
data_shannon_franqueses <- data_shannon_franqueses %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivfranq <- diversity(data_shannon_franqueses)
print(shannondivfranq)## Franqueses-Maig-2021_NS_1 Franqueses-Maig-2021_NS_2 Franqueses-Maig-2021_NS_3
## 1.0789922 0.8675632 0.6931472
## Franqueses-Maig-2021_NS_4 Franqueses-Maig-2021_NS_5 Franqueses-Maig-2021_S_1
## 2.2739657 0.7549968 1.3321790
## Franqueses-Maig-2021_S_2 Franqueses-Maig-2021_S_4
## 1.0986123 1.5607104
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 26 52.32292 16.96108 42.625 8.605049 53.01786 33.13315 4.330973 8
Corba i extrapolació per Sabadell
#Crear una taula de l'abundancia per Mostreig, Tractament i Trampa, amb dades filtrades per Sabadell
abundacia_trampa_sabadell <- dades_finals %>%
filter(Zona == 'Sabadell') %>%
group_by(Mostreig, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament', 'Trampa'. You can
## override using the `.groups` argument.
#Crear una taula de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_auxsaba <- cast(abundacia_trampa_sabadell, Mostreig + Tractament + Trampa ~ ID, value='Numero_total', FUN=mean)
taula_auxsaba <- as.data.frame(taula_auxsaba)
taula_auxsaba[is.na(taula_auxsaba)] <- 0
# Agrupar Mostreig amb Tractament i Trampa
data_shannon_sabadell <- taula_auxsaba %>%
mutate(Agrupacio = paste(Mostreig, Tractament, Trampa, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament, -Trampa)
# Establir tots els valors de la taula com a numèrics: eliminar la columna Agrupacio i establir-la com a identificador de les files
rownames(data_shannon_sabadell) <- data_shannon_sabadell$Agrupacio
data_shannon_sabadell <- data_shannon_sabadell %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# Càlcul de la diversitat de Shannon per cada Mostreig, Tractament i Trampa
shannondivsaba <- diversity(data_shannon_sabadell)
print(shannondivsaba)## Sabadell-Juny-2021_NS_1 Sabadell-Juny-2021_NS_4 Sabadell-Juny-2021_NS_7
## 1.7478681 0.9502705 1.7917595
## Sabadell-Juny-2021_NS_8 Sabadell-Juny-2021_NS_11 Sabadell-Juny-2021_S_2
## 1.4750763 1.3862944 0.6365142
## Sabadell-Juny-2021_S_3 Sabadell-Juny-2021_S_5 Sabadell-Juny-2021_S_6
## 0.5623351 1.4995094 1.0986123
## Sabadell-Juny-2021_S_9 Sabadell-Juny-2021_S_10 Sabadell-Juny-2021_S_12
## 0.6931472 0.5004024 1.2852930
## Species chao chao.se jack1 jack1.se jack2 boot boot.se n
## All 23 45.45833 17.10731 35.83333 5.316379 44.4697 28.44918 2.631169 12
T-student per comparar la diversitat de Shannon entre mostrejos i tractaments
#Convertir la taula de shannondiv en dataframe
shannon_ttest <- as.data.frame(shannondiv)
# Convertir l'ID de les files (mostreig + tractament + número de trampa) en una variable (nova columna)
shannon_ttest$row <- rownames(shannon_ttest)
# De la nova variable, separar en dues columnes el mostreig i eltractament (i no utilitzem el número de trampa)
shannon_ttest <- shannon_ttest %>% mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2], Mostreig = str_split(row, "_", simplify = TRUE)[,1]) %>% dplyr::select(-row)
# Prova t-student i boxplots
shannonT <- shannon_ttest %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(shannondiv[Tractament=='S']),
mean_NS = mean(shannondiv[Tractament=='NS']),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(shannondiv ~ Tractament)$p.value,
NA),.groups = "drop")
shannonT## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 1.02 0.964 0.844
## 2 Campus-Jun-2021 1.74 1.86 0.637
## 3 Campus-juliol-2020 1.08 1.54 0.0840
## 4 Campus-maig-2020 NaN 1.56 NA
## 5 Franqueses-Maig-2021 1.33 1.13 0.566
## 6 Sabadell-Juny-2021 0.897 1.47 0.0232
ggplot(shannon_ttest, aes(x = factor(Mostreig), y = shannondiv, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Diversitat de Shannon x Zona i Tractaments",
x = "Mostreig",
y = "Diversitat Shannon")Només en el mostreig de Sabadell la diversitat de Shannon és significativament diferent (p valor < 0.05): l’índex de Shannon és major al tractament NS.
Diversitat de Shannon d’abelles (per mostreig)
# Preparar les dades: crear taula de dades d'abundància només per abelles
abundanciabelles_ID <- dades_finals %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)
#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux <- cast(abundanciabelles_ID, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux <- as.data.frame(taula_aux)
taula_aux[is.na(taula_aux)] <- 0
#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abelles <- taula_aux %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)
#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abelles) <- data_shannon_abelles$Agrupacio
data_shannon_abelles <- data_shannon_abelles %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
#View(data_shannon_abelles)# Càlcul diversitat de Shannon
shannondiv_abelles <- diversity(data_shannon_abelles)
print(shannondiv_abelles)## Campus-juliol-2020_NS Campus-juliol-2020_S Campus-Juliol-2021_NS
## 2.611229 1.905899 2.090031
## Campus-Juliol-2021_S Campus-Jun-2021_NS Campus-Jun-2021_S
## 1.806531 2.701045 2.080060
## Campus-maig-2020_NS Franqueses-Maig-2021_NS Franqueses-Maig-2021_S
## 2.199104 2.239746 1.560710
## Sabadell-Juny-2021_NS Sabadell-Juny-2021_S
## 2.045357 1.884985
Corbes d’acumulació d’espècies per mostreig
Corba d’acumulació. Riquesa en funció del nº de trampes (specaccum)
# A la taula data_shannon_MTTR extreure una nova columna que indiqui el mostreig
metadata_shannon <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
# View(metadata_shannon)
# Subset each mostreig into its own dataframe (i eliminar la columna de "mostreig" perquè totes les dades de la taula siguin numèriques)
metadata_shannon %>% filter(Zona == 'Campus') %>% dplyr::select(-Zona)-> Campus
metadata_shannon %>% filter(Zona == 'Franqueses') %>% dplyr::select(-Zona) -> Franqueses
metadata_shannon %>% filter(Zona == 'Llig') %>% dplyr::select(-Zona) -> Llig
metadata_shannon %>% filter(Zona == 'Mura') %>% dplyr::select(-Zona) -> Mura
metadata_shannon %>% filter(Zona == 'Sabadell') %>% dplyr::select(-Zona) -> Sabadell
metadata_shannon %>% filter(Zona == 'TV') %>% dplyr::select(-Zona) -> TV7041
# Corba d'acumulació d'espècies per cada mostreig
corba1 = specaccum(Campus, method = "random")
corba2 = specaccum(Franqueses, method = "random")
corba3 = specaccum(Sabadell, method = "random")
# Creating data frames for each set of data
data1 <- data.frame(Sites = corba1$sites, Richness = corba1$richness, SD = corba1$sd)
data1$Zona <- "Campus"
data2 <- data.frame(Sites = corba2$sites, Richness = corba2$richness, SD = corba2$sd)
data2$Zona <- "Franqueses"
data3 <- data.frame(Sites = corba3$sites, Richness = corba3$richness, SD = corba3$sd)
data3$Zona <- "Sabadell"
# Combine all data frames into one
dades_combinades <- rbind(data1, data2, data3)
library(ggplot2)
# Gráfico de dispersión con líneas y error bars
ggplot(dades_combinades, aes(x = Sites, y = Richness, color = Zona)) +
geom_point() +
geom_line(aes(group = Zona), alpha = 0.7) + # Líneas
# geom_ribbon(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD, fill = Mostreig), alpha = 0.01) + # Incertidumbre
geom_errorbar(aes(ymin = Richness - 2 * SD, ymax = Richness + 2 * SD), width = 0.2, alpha = 0.5) + # Barras de error
ylim(0,NA)+
labs(x = "Nº trampes", y = "Riquesa", color = "Zona") +
theme_minimal()Corba de rarefracció. Riquesa en funció del nombre d’individus (rarefy)
dades_agrupades <- metadata_shannon %>%
group_by(Zona) %>%
summarise(across(everything(), sum)) %>%
as.data.frame()
#View (dades_agrupades)
# Buscar el nombre mínim d'observacions en una zona (és 84, a TV juliol)
abundacia_trampa <- dades_finals %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
## # A tibble: 564 × 7
## # Groups: Zona, Mes, Any, Tractament, Trampa [101]
## Zona Mes Any Tractament Trampa ID Numero_total
## <chr> <chr> <dbl> <chr> <dbl> <chr> <int>
## 1 Campus Juliol 2021 NS 2 Lasioglossum glabriusculum 1
## 2 Campus Juliol 2021 NS 4 Curculionidae 1
## 3 Campus Juliol 2021 NS 4 Lasioglossum malachurum 1
## 4 Campus Juliol 2021 NS 4 Pompilidae 1
## 5 Campus Juliol 2021 NS 4 Pyronia cecilia 1
## 6 Campus Juliol 2021 NS 6 Lasioglossum glabriusculum 1
## 7 Campus Juliol 2021 NS 6 Lasioglossum politum 3
## 8 Campus Juliol 2021 NS 6 Pompilidae 1
## 9 Campus Juliol 2021 NS 8 Chrysididae 1
## 10 Campus Juliol 2021 NS 8 Hylaeus sp. 1
## # ℹ 554 more rows
min_n <- abundacia_trampa %>%
group_by(Zona) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dades_agrupades) <- dades_agrupades$Zona
dades_agrupades <- dades_agrupades[,-1]
dades_agrupades <- dades_agrupades %>%
mutate_all(as.numeric)
vegan::rarefy(dades_agrupades,sample=min_n)## Campus Franqueses Sabadell
## 28.12144 26.00000 19.50913
## attr(,"Subsample")
## [1] 59
Corba d’interpolació: riquesa en funció del nombre d’individus (iNEXT). (??)
library(iNEXT)
transposades <- t(dades_agrupades)
D_abund <- iNEXT (transposades, datatype = 'abundance')
plot (D_abund)Corbes de rarefracció d’abelles per cada tractament
# Extreure la columna de tractament a la taula d'abundàncies per ID de cada mostreig
dades_shannon_abelles <- data_shannon_abelles %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])
# Agrupar per tractament
dades_agrupades_tractament <- dades_shannon_abelles %>%
group_by(Tractament) %>%
summarise(across(everything(), sum)) %>%
as.data.frame()
#View (dades_agrupades_tractament)
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abelles <- dades_abelles %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abelles)
min_n <- abundacia_trampa_abelles %>%
group_by(Tractament) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dades_agrupades_tractament) <- dades_agrupades_tractament$Tractament
dades_agrupades_tractament <- dades_agrupades_tractament[,-1]
dades_agrupades_tractament <- dades_agrupades_tractament %>%
mutate_all(as.numeric)
vegan::rarefy(dades_agrupades_tractament,sample=min_n)## NS S
## 59.89088 37.00000
## attr(,"Subsample")
## [1] 321
Corbes de rarefracció d’abelles al Campus per tractament
# Extreure la columna Zona de dades_shannon_abelles i filtrar per Campus
dades_shannon_abellescampus <- dades_shannon_abelles %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
dades_shannon_abellescampus <- dades_shannon_abellescampus[dades_shannon_abellescampus$Zona == "Campus", ]
# Agrupar per tractament
dadescampus_agrupades_tractament <- dades_shannon_abellescampus %>%
group_by(Tractament) %>%
summarise(across(where(is.numeric), sum)) %>%
as.data.frame()
#View (dadescampus_agrupades_tractament)
# Buscar el nombre mínim d'observacions en una zona
abundacia_trampa_abellescampus <- dades_abelles %>%
group_by(Zona, Mes, Any, Tractament, Trampa, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Zona', 'Mes', 'Any', 'Tractament',
## 'Trampa'. You can override using the `.groups` argument.
#View (abundacia_trampa_abellescampus)
min_n <- abundacia_trampa_abellescampus %>%
group_by(Tractament) %>%
summarize (n_especimens = sum(Numero_total)) %>%
summarize(min=min(n_especimens)) %>%
pull(min)
rownames(dadescampus_agrupades_tractament) <- dadescampus_agrupades_tractament$Tractament
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament[,-1]
dadescampus_agrupades_tractament <- dadescampus_agrupades_tractament %>%
mutate_all(as.numeric)
vegan::rarefy(dadescampus_agrupades_tractament,sample=min_n)## Warning in vegan::rarefy(dadescampus_agrupades_tractament, sample = min_n):
## requested 'sample' was larger than smallest site maximum (266)
## NS S
## 54.01949 35.00000
## attr(,"Subsample")
## [1] 321
PROVES T-STUDENT
Comparació de la riquesa per mostrejos i tractaments
Riquesa total
riquesa_x_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Riquesa ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 4.33 3.1 0.303
## 2 Campus-Jun-2021 7.3 9.5 0.190
## 3 Campus-juliol-2020 3.89 5.92 0.0281
## 4 Campus-maig-2020 NaN 6.55 NA
## 5 Franqueses-Maig-2021 4 4.6 0.743
## 6 Sabadell-Juny-2021 3.43 4.8 0.201
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Riquesa per Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")Generalment es compleix un patró de major riquesa al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0,0281)
Riquesa d’abelles
# Filtrar les dades d'abelles
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_abelles %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_abelles <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Riquesa ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_abelles## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 3.89 2.1 0.0847
## 2 Campus-Jun-2021 4.3 5.4 0.338
## 3 Campus-juliol-2020 2.78 4.46 0.00961
## 4 Campus-maig-2020 NaN 3.84 NA
## 5 Franqueses-Maig-2021 2 3 0.341
## 6 Sabadell-Juny-2021 3.14 3.8 0.568
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Riquesa d'abelles per Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")Generalment es compleix un patró de major riquesa d’abelles al tractament NS (excepte a Campus-Juliol-2021). Al mostreig del Campus-Juliol_2020 aquesta diferència és significativa (p-valor= 0.0096)
Riquesa de coleòpters
# Filtrar les dades per coleòpters
dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 29 Campus Juliol 7 2021 4 NS <NA>
## 2 158 Campus Juliol 7 2021 10 NS <NA>
## 3 183 Campus Juliol 7 2021 5 S <NA>
## 4 101 Franqueses Maig 5 2021 2 NS <NA>
## 5 150 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## 6 151 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_coleoptera %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_coleoptera <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Riquesa ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_coleoptera## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 2.88 3.7 0.357
## 3 Campus-juliol-2020 1.33 1.25 0.851
## 4 Campus-maig-2020 NaN 2.44 NA
## 5 Franqueses-Maig-2021 1.5 1 0.500
## 6 Sabadell-Juny-2021 1 1.33 0.423
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Riquesa x Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN COLEÒPTERS.
Riquesa de sírfids
# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 226 Campus Juliol 7 2021 7 S <NA>
## 2 265 Campus Juliol 7 2021 14 NS <NA>
## 3 277 Franqueses Maig 5 2021 4 NS <NA>
## 4 278 Franqueses Maig 5 2021 4 NS <NA>
## 5 279 Franqueses Maig 5 2021 4 NS <NA>
## 6 280 Franqueses Maig 5 2021 4 NS Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
riquesa_x_trampa <- dades_sirfids %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
proves_t_sirfids <- riquesa_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Riquesa[Tractament == "S"]),
mean_NS = mean(Riquesa[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Riquesa ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_sirfids## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <lgl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 1 1 NA
## 3 Campus-juliol-2020 3 1 NA
## 4 Campus-maig-2020 NaN 1.18 NA
## 5 Franqueses-Maig-2021 1 3 NA
## 6 Sabadell-Juny-2021 NaN 1 NA
ggplot(riquesa_x_trampa, aes(x = factor(Mostreig), y = Riquesa, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Riquesa x Zona i Tractaments",
x = "Mostreig",
y = "Riquesa")Hi ha poques dades per obtenir resultats del t-test.
Comparació de l’abundància per mostrejos i tractaments
Abundància total
abundancia_x_trampa <- dades_finals %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_abtotal <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Abundancia ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_abtotal## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 14.1 4.5 0.130
## 2 Campus-Jun-2021 11.6 18.5 0.113
## 3 Campus-juliol-2020 8.89 8.54 0.904
## 4 Campus-maig-2020 NaN 15.4 NA
## 5 Franqueses-Maig-2021 4.67 9 0.178
## 6 Sabadell-Juny-2021 7.29 5.8 0.606
ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")No hi ha cap resultat significatiu (p-valor<0.05). No hi ha un patró clar de com varia l’abundància en funció del tractament.
Abundància d’abelles
# Filtrar les dades d'abelles
dades_abelles <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(dades_abelles) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_abelles %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_ababelles <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = ifelse(length(unique(Tractament)) == 2,
t.test(Abundancia ~ Tractament)$p.value,
NA)
,.groups = "drop")
proves_t_ababelles## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 13.7 3.5 0.108
## 2 Campus-Jun-2021 7.6 8.9 0.589
## 3 Campus-juliol-2020 7.44 7.08 0.896
## 4 Campus-maig-2020 NaN 8.68 NA
## 5 Franqueses-Maig-2021 2 4.4 0.0902
## 6 Sabadell-Juny-2021 7 4.8 0.454
ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")A Franqueses l’abundància d’abelles és significativament superior al tractament NS (p-valor=0.09015). En la resta de mostrejos el patró no és clar.
Abundància de coleòpters
# Filtrar les dades per coleòpters
dades_coleoptera <- dades_finals[dades_finals$Ordre == "Coleoptera", ]
head(dades_coleoptera) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 29 Campus Juliol 7 2021 4 NS <NA>
## 2 158 Campus Juliol 7 2021 10 NS <NA>
## 3 183 Campus Juliol 7 2021 5 S <NA>
## 4 101 Franqueses Maig 5 2021 2 NS <NA>
## 5 150 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## 6 151 Franqueses Maig 5 2021 5 NS Chiqui bupresti…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_coleoptera %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n(),.groups = "drop")
proves_t_abcoleoptera <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Abundancia ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_abcoleoptera## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <dbl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 4.12 9.2 0.175
## 3 Campus-juliol-2020 1.33 1.25 0.851
## 4 Campus-maig-2020 NaN 6.61 NA
## 5 Franqueses-Maig-2021 2 4.67 0.496
## 6 Sabadell-Juny-2021 1 1.33 0.423
ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")
ERROR: NOT ENOUGH OBSERVATIONS. ALGUNS MOSTREJOS NO TENEN
COLEÒPTERS.
Abundància de sírfids
# Filtrar les dades per sírfids
dades_sirfids <-dades_finals %>% filter(Familia=="Syrphidae")
head(dades_sirfids) ## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 226 Campus Juliol 7 2021 7 S <NA>
## 2 265 Campus Juliol 7 2021 14 NS <NA>
## 3 277 Franqueses Maig 5 2021 4 NS <NA>
## 4 278 Franqueses Maig 5 2021 4 NS <NA>
## 5 279 Franqueses Maig 5 2021 4 NS <NA>
## 6 280 Franqueses Maig 5 2021 4 NS Episyrphus balt…
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Proves t-student i plots
abundancia_x_trampa <- dades_sirfids %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Abundancia = n_distinct(ID),.groups = "drop")
proves_t_absirfids <- abundancia_x_trampa %>%
group_by(Mostreig) %>%
summarise(
mean_S = mean(Abundancia[Tractament == "S"]),
mean_NS = mean(Abundancia[Tractament == "NS"]),
t_test_result = tryCatch(
t.test(Abundancia ~ Tractament)$p.value,
error = function(e) NA
)
,.groups = "drop")
proves_t_absirfids## # A tibble: 6 × 4
## Mostreig mean_S mean_NS t_test_result
## <chr> <dbl> <dbl> <lgl>
## 1 Campus-Juliol-2021 1 1 NA
## 2 Campus-Jun-2021 1 1 NA
## 3 Campus-juliol-2020 3 1 NA
## 4 Campus-maig-2020 NaN 1.18 NA
## 5 Franqueses-Maig-2021 1 3 NA
## 6 Sabadell-Juny-2021 NaN 1 NA
ggplot(abundancia_x_trampa, aes(x = factor(Mostreig), y = Abundancia, fill = Tractament, color = Tractament)) +
geom_boxplot(width = 0.7, outlier.shape = NA) +
facet_wrap(~Mostreig, scales = "free") +
labs(title = "Boxplot Abundancia x Zona i Tractaments",
x = "Mostreig",
y = "Abundancia")No hi ha prous dades per obtenir bons resultats de la prova t-student.
DIAGRAMES D’EULER:
Diagrames d’Euler totals per cada zona
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_euler <- as.data.frame.matrix(table(interaction(dades_finals$Zona,dades_finals$ID),dades_finals$Tractament)) %>% mutate(`NS&S`=pmin(NS,S))
head(taula_aux_euler)## NS S NS&S
## Campus.Acmaeodera cylindrica 5 1 1
## Franqueses.Acmaeodera cylindrica 0 0 0
## Sabadell.Acmaeodera cylindrica 0 0 0
## Campus.Amegilla sp. 1 0 0
## Franqueses.Amegilla sp. 0 0 0
## Sabadell.Amegilla sp. 0 0 0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_especies <- taula_aux_euler %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_euler)))
taula_euler_especies <- taula_euler_especies %>% group_by(Zona) %>% summarise( S = list(unique(Especie[ S > 0 ])),
NS = list(unique(Especie[ NS > 0 ])),
`NS&S` = list(unique(Especie[ `NS&S` > 0 ])) )
head(taula_euler_especies)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <list> <list> <list>
## 1 Campus <chr [59]> <chr [98]> <chr [39]>
## 2 Franqueses <chr [11]> <chr [20]> <chr [5]>
## 3 Sabadell <chr [12]> <chr [16]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_especies, file.path(directorio, "2.0_Euler_Especies.xlsx"))# Agrupació per mostrejos (Zona, Mes, Any)
euler_data <- taula_euler_especies %>%
mutate(NS = map_dbl(NS, length),
S = map_dbl(S, length),
`NS&S` = map_dbl(`NS&S`, length))
head(euler_data)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <dbl> <dbl> <dbl>
## 1 Campus 59 98 39
## 2 Franqueses 11 20 5
## 3 Sabadell 12 16 5
# Creació del diagrama d'Euler
apply_euler <- function(row){
euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}
euler_data <- euler_data %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))# Plot del diagrama d'Euler
generate_euler_plots <- function(euler_result,zona){
plot(euler_result,main=paste("Zona: ",zona),quantities = list(type = c("percent", "counts")))}
pmap(list(euler_data$euler_result,euler_data$Zona), generate_euler_plots)## [[1]]
##
## [[2]]
##
## [[3]]
Diagrames d’Euler d’abelles (amb dades_abelles) per cada zona
dades_abelles <- dades_abelles[!(dades_abelles$Mes == "maig" & dades_abelles$Any == 2020 & dades_abelles$Zona == "Campus"), ]
# Preparació de les dades: càlcul de l'abundància de cada espècie (ID) a cada mostreig i tractament (NS, S i intersecció)
taula_aux_eulerabelles <- as.data.frame.matrix(table(interaction(dades_abelles$Zona,dades_abelles$ID),dades_abelles$Tractament)) %>% mutate(`NS&S`=pmin(NS,S))
head(taula_aux_eulerabelles)## NS S NS&S
## Campus.Amegilla sp. 1 0 0
## Franqueses.Amegilla sp. 0 0 0
## Sabadell.Amegilla sp. 0 0 0
## Campus.Andrena cinerea 0 0 0
## Franqueses.Andrena cinerea 1 0 0
## Sabadell.Andrena cinerea 0 0 0
# Reorganitzar la taula auxiliar perquè mostri el nom de les espècies (ID) trobades en cada mostreig i tractament
taula_euler_abelles <- taula_aux_eulerabelles %>% mutate(Especie = sub("^[^.]+\\.", "", rownames(.)), Zona = sub("\\..*","",rownames(taula_aux_eulerabelles)))
taula_euler_abelles <- taula_euler_abelles %>% group_by(Zona) %>% summarise( S = list(unique(Especie[ S > 0 ])),
NS = list(unique(Especie[ NS > 0 ])),
`NS&S` = list(unique(Especie[ `NS&S` > 0 ])) )
head(taula_euler_abelles)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <list> <list> <list>
## 1 Campus <chr [35]> <chr [47]> <chr [21]>
## 2 Franqueses <chr [5]> <chr [12]> <chr [3]>
## 3 Sabadell <chr [10]> <chr [12]> <chr [5]>
# Descarregar les llistes d'espècies en un excel
write.xlsx(taula_euler_abelles, file.path(directorio, "Euler_Abelles.xlsx"))# Agrupació per mostrejos (Zona, Mes, Any)
euler_data_abelles <- taula_euler_abelles %>%
mutate(NS = map_dbl(NS, length),
S = map_dbl(S, length),
`NS&S` = map_dbl(`NS&S`, length))
head(euler_data_abelles)## # A tibble: 3 × 4
## Zona S NS `NS&S`
## <chr> <dbl> <dbl> <dbl>
## 1 Campus 35 47 21
## 2 Franqueses 5 12 3
## 3 Sabadell 10 12 5
# Creació del diagrama d'Euler
apply_euler <- function(row){
euler(c("NS"=row$NS,"S"=row$S,"NS&S"=row$`NS&S`))
}
euler_data_abelles <- euler_data_abelles %>% rowwise() %>% mutate(euler_result= list(apply_euler(cur_data())))# Plot del diagrama d'Euler
generate_euler_plots <- function(euler_result,zona){
plot(euler_result,main=paste("Zona: ",zona),quantities = list(type = c("percent", "counts")))}
pmap(list(euler_data_abelles$euler_result,euler_data_abelles$Zona), generate_euler_plots)## [[1]]
##
## [[2]]
##
## [[3]]
ANÀLISIS DE COMPONENTS PRINCIPALS (PCA)
PCA total
# Preparar dades: taula de dades Shannon per mostreig i tractament (MT). Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig
pca_data <- data_shannon_MT
# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes <- dades_finals %>% group_by(Mostreig, Tractament) %>% summarise(numtrampes=n_distinct(Trampa)) %>% ungroup() %>% mutate(row=paste(Mostreig,Tractament,sep="_")) %>% dplyr::select(-Mostreig,-Tractament)## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data$row <- rownames(pca_data)
pca_data <- pca_data %>% inner_join(n_trampes, by=join_by(row))
# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columns <-ncol(pca_data)
pca_data[, 1:(num_columns - 2)] <- pca_data[, 1:(num_columns - 2)] / pca_data[, num_columns]
# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data <- pca_data %>% dplyr::select(-numtrampes)
pca_data<- pca_data %>%
mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])# Mostrar la PCA
ggbiplot(pca_comps, labels = pca_data$row, groups = pca_data$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))PCA abelles
# Preparar dades: taula de dades Shannon d'abelles. Normalitzar les dades: nombre d'individus entre el nombre de pantraps del mostreig
pca_data_abelles <- data_shannon_abelles
# Calcular el nombre de trampes per cada mostreig i tractament
n_trampes_abelles <- dades_finals %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament) %>%
summarise(numtrampes = n_distinct(Trampa)) %>%
ungroup() %>%
mutate(row = paste(Mostreig, Tractament, sep = "_")) %>%
dplyr::select(-Mostreig, -Tractament)## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
# Unir les dues taules relacionant les files del mateix mostreig i tractament
pca_data_abelles$row <- rownames(pca_data_abelles)
pca_data_abelles <- pca_data_abelles %>% inner_join(n_trampes_abelles, by=join_by(row))
# Dividir les abundàncies de cada ID entre el nombre de trampes de cada mostreig i tractament
num_columnsabe <-ncol(pca_data_abelles)
pca_data_abelles[, 1:(num_columnsabe - 2)] <- pca_data_abelles[, 1:(num_columnsabe - 2)] / pca_data_abelles[, num_columnsabe]
# Eliminar l'última columna de la taula (numtrampes), ja no cal per fer la PCA
pca_data_abelles <- pca_data_abelles %>% dplyr::select(-numtrampes)
pca_data_abelles<- pca_data_abelles %>%
mutate(Tractament = str_split(row, "_", simplify = TRUE)[, 2])
View(pca_data_abelles)# Fer la PCA
pca_abelles <- prcomp(pca_data_abelles %>% dplyr::select(-row,-Tractament))
summary(pca_abelles)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9981 1.1729 0.61490 0.52897 0.48323 0.38196 0.31303
## Proportion of Variance 0.6025 0.2076 0.05706 0.04223 0.03524 0.02202 0.01479
## Cumulative Proportion 0.6025 0.8101 0.86717 0.90940 0.94464 0.96665 0.98144
## PC8 PC9 PC10 PC11
## Standard deviation 0.2482 0.19318 0.15503 2.126e-16
## Proportion of Variance 0.0093 0.00563 0.00363 0.000e+00
## Cumulative Proportion 0.9907 0.99637 1.00000 1.000e+00
# Mostrar la PCA
ggbiplot(pca_abelles, labels = pca_data_abelles$row, groups = pca_data_abelles$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))PCA abelles Campus (per trampes. No sé ben bé què he fet)
# Preparar dades: taula de dades Shannon d'abelles per mostreig, tractament i trampa (dades NO normalitzades)
pca_data_campus <- data_shannon_MTTR %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona == 'Campus') %>%
mutate(Tractament = str_extract(rownames(.), "(?<=_)[^_]+(?=_)"))
#View(pca_data_campus)
# Fer la PCA
pca_campus <- prcomp(pca_data_campus %>% dplyr::select(-Zona,-Tractament))
# Mostrar la PCA
ggbiplot(pca_campus, labels = pca_data_campus$row, groups = pca_data_campus$Tractament,
ellipse = TRUE, var.axes = FALSE) +
theme_minimal() +
theme(text = element_text(size = 20))NMDS, ANÀLISIS DE SIMILITUD (ANOSIM), ADONIS, I ESPÈCIES INDICADORES
NMDS total
# Preparar dades
NMDS_data <- data_shannon_MT
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDS_data %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))
Tractament_metadata## Tractament color
## Campus-juliol-2020_NS NS green
## Campus-juliol-2020_S S orange
## Campus-Juliol-2021_NS NS green
## Campus-Juliol-2021_S S orange
## Campus-Jun-2021_NS NS green
## Campus-Jun-2021_S S orange
## Campus-maig-2020_NS NS green
## Franqueses-Maig-2021_NS NS green
## Franqueses-Maig-2021_S S orange
## Sabadell-Juny-2021_NS NS green
## Sabadell-Juny-2021_S S orange
# NMDS Total
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS <- metaMDS(NMDS_data, k=2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1574807
## Run 1 stress 0.1574807
## ... Procrustes: rmse 1.586753e-06 max resid 1.972763e-06
## ... Similar to previous best
## Run 2 stress 0.1332528
## ... New best solution
## ... Procrustes: rmse 0.2375972 max resid 0.4227959
## Run 3 stress 0.1332526
## ... New best solution
## ... Procrustes: rmse 0.0002214546 max resid 0.0003933833
## ... Similar to previous best
## Run 4 stress 0.1271718
## ... New best solution
## ... Procrustes: rmse 0.07482385 max resid 0.1553409
## Run 5 stress 0.1927377
## Run 6 stress 0.1332526
## Run 7 stress 0.1576285
## Run 8 stress 0.1271718
## ... Procrustes: rmse 6.804317e-05 max resid 0.0001405358
## ... Similar to previous best
## Run 9 stress 0.3407341
## Run 10 stress 0.1271718
## ... Procrustes: rmse 1.536335e-05 max resid 3.205959e-05
## ... Similar to previous best
## Run 11 stress 0.1332526
## Run 12 stress 0.1576277
## Run 13 stress 0.1271718
## ... Procrustes: rmse 4.444601e-06 max resid 9.072617e-06
## ... Similar to previous best
## Run 14 stress 0.1647707
## Run 15 stress 0.1271718
## ... Procrustes: rmse 1.498339e-05 max resid 3.060948e-05
## ... Similar to previous best
## Run 16 stress 0.1574807
## Run 17 stress 0.1332526
## Run 18 stress 0.1576295
## Run 19 stress 0.1332525
## Run 20 stress 0.1927396
## *** Best solution repeated 4 times
# Plot NMDS total
ordiplot(NMDS,type="n")
ordihull(NMDS,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDS,display="species",col="black",air=0.01)
orditorp(NMDS,display="sites",col=Tractament_metadata$color,
air=0.01,cex=1.25)# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("Tot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()## quartz_off_screen
## 2
Anàlisi de similitud (ANOSIM) Total
##
## Call:
## anosim(x = NMDS_data, grouping = Tractament_metadata$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.06667
## Significance: 0.721
##
## Permutation: free
## Number of permutations: 9999
No signiticatiu (p valor > 0.05): la composició de les comunitats no és significativament diferent entre els tractaments S i NS. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.
ADONIS total
# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata_TOTAL <- NMDS_data %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)
Zona_Tractament_metadata_TOTAL## Zona Tractament
## Campus-juliol-2020_NS Campus NS
## Campus-juliol-2020_S Campus S
## Campus-Juliol-2021_NS Campus NS
## Campus-Juliol-2021_S Campus S
## Campus-Jun-2021_NS Campus NS
## Campus-Jun-2021_S Campus S
## Campus-maig-2020_NS Campus NS
## Franqueses-Maig-2021_NS Franqueses NS
## Franqueses-Maig-2021_S Franqueses S
## Sabadell-Juny-2021_NS Sabadell NS
## Sabadell-Juny-2021_S Sabadell S
NMDStot.dist <- vegdist(wisconsin(sqrt(NMDS_data)), k=2)
adonis(formula = NMDStot.dist~Tractament/Zona, data = Zona_Tractament_metadata_TOTAL, permutations = 1000)$aov.tab## 'adonis' will be deprecated: use 'adonis2' instead
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.3611 0.36111 1.1373 0.10191 0.25874
## Tractament:Zona 4 1.5947 0.39868 1.2556 0.45004 0.03996 *
## Residuals 5 1.5877 0.31753 0.44805
## Total 10 3.5435 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No hi ha diferències significatives entre tractaments (p-valor = 0.25874), però la interacció entre zona i tractament explica un 45% de la variabilitat, i és un resultat significatiu (p-valor=0.03996).
Espècies indicadores dels tractaments S i NS total
indicadores = multipatt(NMDS_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 137
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No hi ha espècies indicadores significatives de cap dels dos tractaments.
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.2
##
## Total number of species: 137
## Selected number of species: 6
## Number of species associated to 1 group: 6
##
## List of species associated to each combination:
##
## Group NS #sps. 3
## stat p.value
## Anthaxia sp. 0.422 0.182
## Eucera elongatula 0.422 0.181
## Ceratina parvula 0.391 0.178
##
## Group S #sps. 3
## stat p.value
## Vespula germanica 0.500 0.186
## Lasioglossum glabriusculum 0.430 0.184
## Lasioglossum griseolum 0.385 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Anthaxia sp., Eucera elongatula, Ceratina parvula), i 3 espècies indicadores de S (Vespula germanica, Lasioglossum gabriusculum i Lasioglossum griseolum). Però realment no son significatives.
NMDS Campus
# Filtrar les dades del Campus i preparar les dades
NMDS_data_campus <- NMDS_data %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona=="Campus") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDS_data_campus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))
# NMDS_data_campus# NMDS Campus
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDScampus <- metaMDS(NMDS_data_campus, k=2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.03999251
## Run 1 stress 0.03999252
## ... Procrustes: rmse 3.432746e-05 max resid 4.816299e-05
## ... Similar to previous best
## Run 2 stress 0.03999251
## ... New best solution
## ... Procrustes: rmse 1.138143e-05 max resid 1.575934e-05
## ... Similar to previous best
## Run 3 stress 0.03999257
## ... Procrustes: rmse 0.0001319502 max resid 0.0001794008
## ... Similar to previous best
## Run 4 stress 0.2570537
## Run 5 stress 0.03999251
## ... Procrustes: rmse 1.730368e-05 max resid 2.266611e-05
## ... Similar to previous best
## Run 6 stress 0.04217697
## Run 7 stress 0.042177
## Run 8 stress 0.2570032
## Run 9 stress 0.04581483
## Run 10 stress 0.04581482
## Run 11 stress 0.2570537
## Run 12 stress 0.04318443
## Run 13 stress 0.2570537
## Run 14 stress 0.04318432
## Run 15 stress 0.0432258
## Run 16 stress 0.03999251
## ... Procrustes: rmse 1.482389e-05 max resid 2.245167e-05
## ... Similar to previous best
## Run 17 stress 0.293306
## Run 18 stress 0.04581486
## Run 19 stress 0.04625739
## Run 20 stress 0.1960094
## *** Best solution repeated 4 times
# Plot NMDS Campus
ordiplot(NMDScampus,type="n")
ordihull(NMDScampus,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDScampus,display="species",col="black",air=0.01)
orditorp(NMDScampus,display="sites",col=Tractament_metadata$color,
air=0.01,cex=1.25)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("output_plot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()Anàlisi de similitud (ANOSIM) Campus
ano = anosim(NMDS_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
##
## Call:
## anosim(x = NMDS_data_campus, grouping = Tractament_metadata$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.03704
## Significance: 0.42857
##
## Permutation: free
## Number of permutations: 5039
No signiticatiu (p valor = 0.42857): la composició de les comunitats de S i NS al Campus no és significativament diferent entre els dos tractaments (S i NS).
ADONIS Campus
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata_CAMPUS <- NMDS_data_campus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
Tractament_metadata_CAMPUS## Tractament
## Campus-juliol-2020_NS NS
## Campus-juliol-2020_S S
## Campus-Juliol-2021_NS NS
## Campus-Juliol-2021_S S
## Campus-Jun-2021_NS NS
## Campus-Jun-2021_S S
## Campus-maig-2020_NS NS
NMDScamp.dist <- vegdist(wisconsin(sqrt(NMDS_data_campus)), k=2)
adonis(formula = NMDScamp.dist~Tractament, data = Tractament_metadata_CAMPUS, permutations = 1000)$aov.tab## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.27747 0.27747 0.87044 0.14828 0.5924
## Residuals 5 1.59387 0.31877 0.85172
## Total 6 1.87134 1.00000
El tractament només explica un 14,8% de la variabilitat, i el resultat no és significatiu.
Espècies indicadores dels tractaments S i NS Campus
indicadores = multipatt(NMDS_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 137
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No hi ha espècies indicadores significatives de cap dels dos tractaments.
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.2
##
## Total number of species: 137
## Selected number of species: 5
## Number of species associated to 1 group: 5
##
## List of species associated to each combination:
##
## Group NS #sps. 3
## stat p.value
## Halictus scabiosae 0.775 0.149
## Anthaxia sp. 0.542 0.147
## Ceratina parvula 0.498 0.147
##
## Group S #sps. 2
## stat p.value
## Lasioglossum glabriusculum 0.731 0.118
## Vespula germanica 0.707 0.145
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Si augmentem el nivell de significació (alfa) a 0,2, hi ha 3 espècies indicadores de NS (Halictus scabiosae, Anthaxia sp., i Ceratina parvula), i 2 espècies indicadores de S (Lasioglossum glabriusculum i Vespula germanica). Però no son resultats significatius.
NMDS d’abelles totals (amb data_shannon_abelles)
# Crear metadades del tractament
NMDSabelles_data <- data_shannon_abelles
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDSabelles_data %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))
Tractament_metadata## Tractament color
## Campus-juliol-2020_NS NS green
## Campus-juliol-2020_S S orange
## Campus-Juliol-2021_NS NS green
## Campus-Juliol-2021_S S orange
## Campus-Jun-2021_NS NS green
## Campus-Jun-2021_S S orange
## Campus-maig-2020_NS NS green
## Franqueses-Maig-2021_NS NS green
## Franqueses-Maig-2021_S S orange
## Sabadell-Juny-2021_NS NS green
## Sabadell-Juny-2021_S S orange
# NMDS abelles
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDSabe <- metaMDS(NMDSabelles_data, k=2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1466713
## Run 1 stress 0.1762197
## Run 2 stress 0.1764663
## Run 3 stress 0.1430526
## ... New best solution
## ... Procrustes: rmse 0.1296798 max resid 0.3307691
## Run 4 stress 0.1660287
## Run 5 stress 0.1350214
## ... New best solution
## ... Procrustes: rmse 0.08009816 max resid 0.2022004
## Run 6 stress 0.1430526
## Run 7 stress 0.1660287
## Run 8 stress 0.1350213
## ... New best solution
## ... Procrustes: rmse 0.000511821 max resid 0.0009819792
## ... Similar to previous best
## Run 9 stress 0.3407333
## Run 10 stress 0.1350213
## ... Procrustes: rmse 0.0003994713 max resid 0.0007627691
## ... Similar to previous best
## Run 11 stress 0.1466713
## Run 12 stress 0.1350216
## ... Procrustes: rmse 0.0005481505 max resid 0.001090682
## ... Similar to previous best
## Run 13 stress 0.1430527
## Run 14 stress 0.1487731
## Run 15 stress 0.1430526
## Run 16 stress 0.1660287
## Run 17 stress 0.1487731
## Run 18 stress 0.211699
## Run 19 stress 0.1466713
## Run 20 stress 0.1350213
## ... Procrustes: rmse 0.0002187493 max resid 0.000402306
## ... Similar to previous best
## *** Best solution repeated 4 times
##
## Call:
## metaMDS(comm = NMDSabelles_data, k = 2)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(NMDSabelles_data))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1350213
## Stress type 1, weak ties
## Best solution was repeated 4 times in 20 tries
## The best solution was from try 8 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(NMDSabelles_data))'
# Plot NMDS abelles
ordiplot(NMDSabe,type="n")
ordihull(NMDSabe,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDSabe,display="species",col="black",air=0.01)
orditorp(NMDSabe,display="sites",col=Tractament_metadata$color,
air=0.01,cex=1.25)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("AbellesTot.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()Anàlisi de similitud (ANOSIM) abelles totals
ano = anosim(NMDSabelles_data, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)
ano##
## Call:
## anosim(x = NMDSabelles_data, grouping = Tractament_metadata$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.04267
## Significance: 0.6117
##
## Permutation: free
## Number of permutations: 9999
No signiticatiu (p valor > 0.6037): la composició de les comunitats és pràcticament idèntica (valor de l’estadístic R proper a 0, i fins i tot negatiu). A més, el resultat no és significatiu. També s’observa al plot de NMDS: els polígons que corresponen als dos tractaments estan superposats. Això probablement sigui per haver fet la prova amb totes les zones de mostreig alhora.
ADONIS abelles totals
# Extreure de la taula el tractament de cada mostreig: una nova columna
Zona_Tractament_metadata <- NMDSabelles_data %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Zona, Tractament)
Zona_Tractament_metadata## Zona Tractament
## Campus-juliol-2020_NS Campus NS
## Campus-juliol-2020_S Campus S
## Campus-Juliol-2021_NS Campus NS
## Campus-Juliol-2021_S Campus S
## Campus-Jun-2021_NS Campus NS
## Campus-Jun-2021_S Campus S
## Campus-maig-2020_NS Campus NS
## Franqueses-Maig-2021_NS Franqueses NS
## Franqueses-Maig-2021_S Franqueses S
## Sabadell-Juny-2021_NS Sabadell NS
## Sabadell-Juny-2021_S Sabadell S
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data)), k=2)
adonis(formula = NMDSabe.dist~Tractament/Zona, data = Zona_Tractament_metadata, permutations = 1000)$aov.tab## 'adonis' will be deprecated: use 'adonis2' instead
## Permutation: free
## Number of permutations: 1000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.3748 0.37477 1.3134 0.11721 0.1299
## Tractament:Zona 4 1.3958 0.34896 1.2229 0.43656 0.0979 .
## Residuals 5 1.4267 0.28534 0.44622
## Total 10 3.1973 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El tractament només explica un 11,72% de la variabilitat observada entre les poblacions de S i NS, i no és un resultat significatiu. La interacció entre el tractament i la zona explica un 43,65% de la variabilitat observada, i és un resultat marginalment significatiu (0,0989, < 0.1).
Espècies d’abelles indicadores dels tractaments S i NS
indicadores = multipatt(NMDSabelles_data, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 80
## Selected number of species: 0
## Number of species associated to 1 group: 0
##
## List of species associated to each combination:
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No hi ha espècies indicadores significatives.
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.2
##
## Total number of species: 80
## Selected number of species: 3
## Number of species associated to 1 group: 3
##
## List of species associated to each combination:
##
## Group NS #sps. 1
## stat p.value
## Eucera elongatula 0.422 0.181
##
## Group S #sps. 2
## stat p.value
## Lasioglossum glabriusculum 0.430 0.184
## Lasioglossum griseolum 0.385 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Si augmentem el nivell de significació (alfa) a 0.2, hi ha 1 espècie indicadora a NS (Eucera elongatula) i 2 espècies indicadores a S (L. glabriusculum i L. griseolum). Però no son resultats significatius.
NMDS abelles al Campus:
# Filtrar les dades del Campus i preparar les dades
NMDSabelles_data_campus <- NMDSabelles_data %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
filter(Zona=="Campus") %>% dplyr::select(-Zona)
# Extreure de la taula el tractament de cada mostreig: una nova columna
Tractament_metadata <- NMDSabelles_data_campus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
# Crear una nova columna a la taula, que associa un color a cada tractament (NS verd, i S blau), per fer el gràfic
Tractament_metadata <- Tractament_metadata %>% mutate(color=ifelse(Tractament=="NS","green","orange"))# NMDS abelles Campus
# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDSabe_Campus <- metaMDS(NMDSabelles_data_campus, k=2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.05473171
## Run 1 stress 0.05473171
## ... Procrustes: rmse 1.536864e-06 max resid 2.917548e-06
## ... Similar to previous best
## Run 2 stress 0.09163904
## Run 3 stress 0.05587716
## Run 4 stress 0.2601494
## Run 5 stress 0.05473171
## ... Procrustes: rmse 3.481031e-06 max resid 6.748059e-06
## ... Similar to previous best
## Run 6 stress 0.05473171
## ... New best solution
## ... Procrustes: rmse 1.075571e-06 max resid 1.942485e-06
## ... Similar to previous best
## Run 7 stress 0.2143136
## Run 8 stress 0.05473171
## ... Procrustes: rmse 3.503063e-06 max resid 6.969015e-06
## ... Similar to previous best
## Run 9 stress 0.05587716
## Run 10 stress 0.1001423
## Run 11 stress 0.05587716
## Run 12 stress 0.05473171
## ... Procrustes: rmse 3.712786e-06 max resid 7.351559e-06
## ... Similar to previous best
## Run 13 stress 0.05587716
## Run 14 stress 0.05473171
## ... Procrustes: rmse 3.16519e-06 max resid 6.1088e-06
## ... Similar to previous best
## Run 15 stress 0.2294233
## Run 16 stress 0.1001423
## Run 17 stress 0.293306
## Run 18 stress 0.05587716
## Run 19 stress 0.1001423
## Run 20 stress 0.1964881
## *** Best solution repeated 4 times
# Plot NMDS
ordiplot(NMDSabe_Campus,type="n")
ordihull(NMDSabe_Campus,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDSabe_Campus,display="species",col="black",air=0.01)
orditorp(NMDSabe_Campus,display="sites",col=Tractament_metadata$color,
air=0.01,cex=1.25)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("AbellesCampus.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()Anàlisi de similitud (ANOSIM) abelles Campus
ano = anosim(NMDSabelles_data_campus, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
##
## Call:
## anosim(x = NMDSabelles_data_campus, grouping = Tractament_metadata$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.03704
## Significance: 0.34286
##
## Permutation: free
## Number of permutations: 5039
La composició de les comunitats entre els dos tractaments és molt similar (valor de l’estadístic d’R molt proper a 0), i el resultat no és significatiu.
ADONIS abelles Campus
# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS1_Tractament_metadata <- NMDSabelles_data_campus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
ADONIS1_Tractament_metadata## Tractament
## Campus-juliol-2020_NS NS
## Campus-juliol-2020_S S
## Campus-Juliol-2021_NS NS
## Campus-Juliol-2021_S S
## Campus-Jun-2021_NS NS
## Campus-Jun-2021_S S
## Campus-maig-2020_NS NS
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDSabelles_data_campus)), k=2)
adonis(formula = NMDSabe.dist~Tractament, data = ADONIS1_Tractament_metadata, permutations = 1000)$aov.tab## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.26326 0.26326 0.9175 0.15505 0.4905
## Residuals 5 1.43466 0.28693 0.84495
## Total 6 1.69792 1.00000
El factor Tractament explica un 15.5% de la variabilitat d’abelles al Campus entre els tractaments S i NS, però aquesta explicació no és significativa.
Espècies d’abelles indicadores dels tractaments S i NS al Campus
indicadores = multipatt(NMDSabelles_data_campus, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.2
##
## Total number of species: 80
## Selected number of species: 2
## Number of species associated to 1 group: 2
##
## List of species associated to each combination:
##
## Group NS #sps. 1
## stat p.value
## Halictus scabiosae 0.775 0.149
##
## Group S #sps. 1
## stat p.value
## Lasioglossum glabriusculum 0.731 0.118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No hi ha espècies indicadores significatives. Si augmenta el nivell de significació alfa = 0.2, sí que apareixen Halictus scabiosae com a indicador del tractament NO segat, i Lasioglossum glabriusculum com a indicador del tractament segat.
NMDS amb presència-absència (no abundàncies) ABELLES CAMPUS
NMDS_data_abellescampus_presabs<-NMDSabelles_data_campus
NMDS_data_abellescampus_presabs[NMDS_data_abellescampus_presabs>0]=1# set seed: si hi ha algun factor d'aleatorietat, fixar una "llavor" perquè sempre surti el mateix resultat (el 7 perquè m'agrada)
set.seed(7)
# k=2: nombre de dimensions
NMDS_presabs <- metaMDS(NMDS_data_abellescampus_presabs, k=2)
ordiplot(NMDS_presabs,type="n")
ordihull(NMDS_presabs,groups=Tractament_metadata$Tractament,draw="polygon",col="grey90",label=F)
orditorp(NMDS_presabs,display="species",col="black",air=0.01)
orditorp(NMDS_presabs,display="sites",col=Tractament_metadata$color,
air=0.01,cex=1.25)
# Save the last plot
g <- recordPlot()
# Replay and save as PNG
width <- 17 # Width in inches
height <- 14 # Height in inches
res <- 500 # Resolution in dpi
# Replay and save as PNG with higher resolution and larger size
png("Presenciabsencia.png", width = width, height = height, units = "in", res = res)
replayPlot(g)
dev.off()Anàlisi de similitud ANOSIM presència-absència
ano = anosim(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, distance = "bray", permutations = 9999)## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
##
## Call:
## anosim(x = NMDS_data_abellescampus_presabs, grouping = Tractament_metadata$Tractament, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: -0.1019
## Significance: 0.62857
##
## Permutation: free
## Number of permutations: 5039
Les composicions de les comunitats dels dos tractaments son molt similars, i el resultat no és significatiu.
ADONIS presència-absència
# Extreure de la taula el tractament de cada mostreig: una nova columna
ADONIS_Tractament_metadata <- NMDS_data_abellescampus_presabs %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>% dplyr::select(Tractament)
ADONIS_Tractament_metadata## Tractament
## Campus-juliol-2020_NS NS
## Campus-juliol-2020_S S
## Campus-Juliol-2021_NS NS
## Campus-Juliol-2021_S S
## Campus-Jun-2021_NS NS
## Campus-Jun-2021_S S
## Campus-maig-2020_NS NS
NMDSabe.dist <- vegdist(wisconsin(sqrt(NMDS_data_abellescampus_presabs)), k=2)
adonis(formula = NMDSabe.dist~Tractament, data = ADONIS_Tractament_metadata, permutations = 1000)$aov.tab## 'adonis' will be deprecated: use 'adonis2' instead
## Set of permutations < 'minperm'. Generating entire set.
## Permutation: free
## Number of permutations: 5039
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tractament 1 0.23205 0.23205 0.97825 0.16363 0.3716
## Residuals 5 1.18603 0.23721 0.83637
## Total 6 1.41808 1.00000
El tractament explica un 16% de la variabilitat entre les comunitats, però el resultat no és significatiu.
Espècies indicadores presència-absència
indicadores = multipatt(NMDS_data_abellescampus_presabs, Tractament_metadata$Tractament, func="r.g", control = how(nperm = 9999))
summary(indicadores, alpha = 0.2)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.2
##
## Total number of species: 80
## Selected number of species: 1
## Number of species associated to 1 group: 1
##
## List of species associated to each combination:
##
## Group NS #sps. 1
## stat p.value
## Halictus scabiosae 0.775 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Només Halictus scabiosae com a indicadora del tractament NS, però el resultat no és significatiu (nivell de significació alfa = 0.2)
MODELS LINEALS MIXTOS D’ANÀLISIS DESCRIPTIVES:
Abundància de Lasioglossums
# Crear matriu amb només les columnes de Lasioglossums
noms_columnes <- colnames(metadata_shannon)
genere <- "Lasioglossum"
columnes_lasioglossum <- noms_columnes[grep(paste0("^", genere), noms_columnes)]
dades_lasioglossum <- metadata_shannon[, columnes_lasioglossum]
# Agrupar totes les columnes en una sola (abundància total de Lasioglossums en cada trampa)
dades_lasioglossum$Abundancia_Total <- rowSums(dades_lasioglossum)
dades_lasioglossum_ok <- data.frame(Abundancia_Total = dades_lasioglossum$Abundancia_Total)
rownames(dades_lasioglossum_ok) <- rownames(dades_lasioglossum)
# View(dades_lasioglossum_ok)
# Extreure columnes de zona i tractament
lasioglossums_zona_tractament <- dades_lasioglossum_ok %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_extract(rownames(.), "(?<=_)[A-Z]+(?=_\\d)"))
# View(lasioglossums_zona_tractament)Model lineal mixt (tractament com a variable explicativa, i zona com a variable aleatòria)
# install.packages("lme4")
# install.packages("Matrix")
model_lasioglossum <- lmer(Abundancia_Total ~ Tractament + (1|Zona), data = lasioglossums_zona_tractament)
summary(model_lasioglossum)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Abundancia_Total ~ Tractament + (1 | Zona)
## Data: lasioglossums_zona_tractament
##
## REML criterion at convergence: 642.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1440 -0.4778 -0.1408 0.1962 7.2804
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 1.557 1.248
## Residual 35.226 5.935
## Number of obs: 101, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.926 1.172 2.051 1.643 0.23908
## TractamentS 3.954 1.226 98.906 3.224 0.00171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.429
Quan el tractament és NS, l’abundància de Lasioglossums és superior a 0 (Intercept: marginalment significatiu). Quan el tractament és S, s’espera que l’abundància de lasioglossums augmenti 3,954 respecte el tractament NS, i és un resultat molt significatiu (p-valor = 0.00171). Els Lasioglossums son més abundants a les zones segades (concorda amb algunes de les espècies “indicadores” de S!)
Model lineal (Tractament com a variable explicativa, sense considerar la zona)
# Convertir la variable Tractament a factor
lasioglossums_zona_tractament$Tractament <- factor(lasioglossums_zona_tractament$Tractament)
# Ajustar el model lineal
model <- lm(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
summary(model)##
## Call:
## lm(formula = Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.421 -2.619 -1.421 1.381 43.579
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.619 0.754 3.473 0.000764 ***
## TractamentS 3.802 1.229 3.093 0.002576 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.985 on 99 degrees of freedom
## Multiple R-squared: 0.08811, Adjusted R-squared: 0.0789
## F-statistic: 9.566 on 1 and 99 DF, p-value: 0.002576
# Crear el gráfico
plot(Abundancia_Total ~ Tractament, data = lasioglossums_zona_tractament, pch = 16)Diversitat de Shannon d’abelles (zona com a variable aleatòria)
# Crear un dataframe de la diversitat de shannon d'abelles (calculada abans)
taula_shannon_abelles <- data.frame(Shannon = shannondiv_abelles)
rownames(taula_shannon_abelles) <- rownames(data_shannon_abelles)
# Extreure dues noves columnes: tractament i zona
taula_shannon_abelles <- taula_shannon_abelles %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1]) %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2])
print(taula_shannon_abelles)## Shannon Zona Tractament
## Campus-juliol-2020_NS 2.611229 Campus NS
## Campus-juliol-2020_S 1.905899 Campus S
## Campus-Juliol-2021_NS 2.090031 Campus NS
## Campus-Juliol-2021_S 1.806531 Campus S
## Campus-Jun-2021_NS 2.701045 Campus NS
## Campus-Jun-2021_S 2.080060 Campus S
## Campus-maig-2020_NS 2.199104 Campus NS
## Franqueses-Maig-2021_NS 2.239746 Franqueses NS
## Franqueses-Maig-2021_S 1.560710 Franqueses S
## Sabadell-Juny-2021_NS 2.045357 Sabadell NS
## Sabadell-Juny-2021_S 1.884985 Sabadell S
# Fer el model
model_shannonabelles <- lmer(Shannon ~ Tractament + (1|Zona), data = taula_shannon_abelles)
summary(model_shannonabelles)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Zona)
## Data: taula_shannon_abelles
##
## REML criterion at convergence: 3.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.12461 -0.79509 -0.01592 0.61598 1.56288
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 0.009103 0.09541
## Residual 0.051690 0.22736
## Number of obs: 11, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2831 0.1121 3.4370 20.373 0.000105 ***
## TractamentS -0.4605 0.1378 7.5671 -3.342 0.011055 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.570
Quan el tractament és NS, el valor de l’índex de Shannon és significativament superior a 0 (Intercept). Quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,4605 respecte el tractament NS. És un resultat significatiu (p-valor = 0,011055). Per tant, la diversitat de Shannon en abelles és significativament inferior en el tractament segat.
# Plot del model
plot(fitted(model_shannonabelles), resid(model_shannonabelles),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Riquesa total
# Preparar les dades
riquesa_aux <- as.data.frame(riquesa)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)
# Extreure columna de la zona
taula_riquesa <- riquesa_aux %>%
mutate(Zona = str_split(rownames(.), "-", simplify = TRUE)[, 1])
print(taula_riquesa)## Tractament riquesa Zona
## Campus-Juliol-2021_NS NS 20 Campus
## Campus-Juliol-2021_S S 21 Campus
## Campus-Jun-2021_NS NS 45 Campus
## Campus-Jun-2021_S S 34 Campus
## Campus-juliol-2020_NS NS 39 Campus
## Campus-juliol-2020_S S 24 Campus
## Campus-maig-2020_NS NS 45 Campus
## Franqueses-Maig-2021_NS NS 20 Franqueses
## Franqueses-Maig-2021_S S 11 Franqueses
## Sabadell-Juny-2021_NS NS 16 Sabadell
## Sabadell-Juny-2021_S S 12 Sabadell
# Fer el model
model_riquesa <- lmer(riquesa ~ Tractament + (1|Zona), data = taula_riquesa)
summary(model_riquesa)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Zona)
## Data: taula_riquesa
##
## REML criterion at convergence: 70.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7954 -0.4064 -0.1774 0.6667 1.1012
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 86.50 9.301
## Residual 74.49 8.631
## Number of obs: 11, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 26.171 6.627 2.885 3.949 0.0311 *
## TractamentS -9.501 5.237 7.267 -1.814 0.1109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.373
TractamentS: quan el tractament és S la riquesa disminueix (-9.501), però no és un resultat significatiu.
# Plot del model
plot(fitted(model_riquesa), resid(model_riquesa),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Riquesa PER TRAMPA d’abelles al Campus (data de mostreig com a factor aleatori)
# Preparar les dades
dades_abellescampus <- dades_abelles[dades_abelles$Zona == "Campus", ]
riquesatrampa_abellescampus <- dades_abellescampus %>%
group_by(Mostreig, Tractament, Trampa) %>%
summarise(Riquesa = n_distinct(ID),.groups = "drop")
riquesa_abellescampus <- as.data.frame(riquesatrampa_abellescampus)
rownames(riquesa_abellescampus) <- paste(riquesa_abellescampus$Mostreig, riquesa_abellescampus$Tractament, riquesa_abellescampus$Trampa, sep = "_")
riquesa_abellescampus <- subset(riquesa_abellescampus, select = -Mostreig)
# Extreure columna de la data de mostreig:
riquesa_abellescampus <- riquesa_abellescampus %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
head(riquesa_abellescampus)## Tractament Trampa Riquesa Data_Mostreig
## Campus-Juliol-2021_NS_2 NS 2 1 Juliol-2021
## Campus-Juliol-2021_NS_4 NS 4 1 Juliol-2021
## Campus-Juliol-2021_NS_6 NS 6 2 Juliol-2021
## Campus-Juliol-2021_NS_8 NS 8 3 Juliol-2021
## Campus-Juliol-2021_NS_10 NS 10 1 Juliol-2021
## Campus-Juliol-2021_NS_12 NS 12 3 Juliol-2021
# Fer el model
model_riquesa_abellescampus <- lmer(Riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_abellescampus)
summary(model_riquesa_abellescampus)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Riquesa ~ Tractament + (1 | Data_Mostreig)
## Data: riquesa_abellescampus
##
## REML criterion at convergence: 265
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8026 -0.7445 -0.1710 0.4959 2.9434
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.6885 0.8297
## Residual 4.4396 2.1070
## Number of obs: 61, groups: Data_Mostreig, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.0378 0.6041 2.8500 6.684 0.00801 **
## TractamentS -0.3864 0.5427 57.1316 -0.712 0.47932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.414
En el tractament S la riquesa d’abelles per trampa disminueix, però no és significatiu.
# Plot del model
plot(fitted(model_riquesa_abellescampus), resid(model_riquesa_abellescampus),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Riquesa PER MOSTREIG d’abelles al Campus (data de mostreig com a factor aleatori)
# Preparar les dades
riquesa_aux <- dades_abellescampus %>%
group_by(Mostreig, Tractament) %>%
summarise(riquesa = n_distinct(ID))## `summarise()` has grouped output by 'Mostreig'. You can override using the
## `.groups` argument.
riquesa_aux <- as.data.frame(riquesa_aux)
rownames(riquesa_aux) <- paste(riquesa_aux$Mostreig, riquesa_aux$Tractament, sep = "_")
riquesa_aux <- subset(riquesa_aux, select = -Mostreig)
# Extreure columna de la data de mostreig:
riquesa_aux <- riquesa_aux %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
riquesa_aux## Tractament riquesa Data_Mostreig
## Campus-Juliol-2021_NS NS 12 Juliol-2021
## Campus-Juliol-2021_S S 17 Juliol-2021
## Campus-Jun-2021_NS NS 25 Jun-2021
## Campus-Jun-2021_S S 19 Jun-2021
## Campus-juliol-2020_NS NS 27 juliol-2020
## Campus-juliol-2020_S S 14 juliol-2020
# Fer el model
model_riquesa_abellescampus2 <- lmer(riquesa ~ Tractament + (1|Data_Mostreig), data = riquesa_aux)## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: riquesa ~ Tractament + (1 | Data_Mostreig)
## Data: riquesa_aux
##
## REML criterion at convergence: 27.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5484 -0.3180 0.2212 0.5530 0.9401
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.00 0.000
## Residual 36.33 6.028
## Number of obs: 6, groups: Data_Mostreig, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.333 3.480 4.000 6.130 0.00359 **
## TractamentS -4.667 4.922 4.000 -0.948 0.39672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.707
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
TractamentS: quan el tractament és S la riquesa disminueix (en 5.583), però NO és un resultat significatiu (p-valor = 0.24677).
# Plot del model
plot(fitted(model_riquesa_abellescampus2), resid(model_riquesa_abellescampus2),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Diversitat (Shannon) d’abelles al Campus (data de mostreig com a factor aleatori)
# Crear taula auxiliar: abundàncies d'abelles per ID, per cada mostreig del Campus
# Filtrar dades per Campus
dades_abellescampus <- dades_finals[dades_finals$Zona == "Campus", ]
# Crear taula de dades d'abundància només per abelles
abundanciabelles_ID2 <- dades_abellescampus %>%
mutate(Mostreig = paste(Zona, Mes, Any, sep = "-")) %>%
filter(Ordre == 'HymenopteraAbella') %>%
group_by(Mostreig, Tractament, ID) %>%
summarise(Numero_total = n())## `summarise()` has grouped output by 'Mostreig', 'Tractament'. You can override
## using the `.groups` argument.
# View(abundanciabelles_ID)
#Crear una taula auxiliar de l'abundància de cada ID (columnes) segons el Mostreig i el tractament (files). És a dir, reorganitzar la taula d'abundancia_ID.
taula_aux2 <- cast(abundanciabelles_ID2, Mostreig + Tractament ~ ID, value='Numero_total', FUN=mean)
taula_aux2 <- as.data.frame(taula_aux2)
taula_aux2[is.na(taula_aux2)] <- 0
#Agrupar Mostreig i tractament en una nova variable (Agrupacio)
data_shannon_abellescampus <- taula_aux2 %>%
mutate(Agrupacio = paste(Mostreig, Tractament, sep = "_")) %>%
ungroup() %>%
dplyr::select(-Mostreig, -Tractament)
#Definir els rownames amb la variable Agrupacio (mostreig i tractament)
rownames(data_shannon_abellescampus) <- data_shannon_abellescampus$Agrupacio
data_shannon_abellescampus <- data_shannon_abellescampus %>% dplyr::select(-Agrupacio) %>% mutate_all(as.numeric)
# head(data_shannon_abellescampus)# Càlcul de l'índex de Shannon i preparació del dataframe
shannondiv_abellescampus <- diversity(data_shannon_abellescampus)
# Crear un dataframe de la diversitat de shannon d'abelles
taula_shannon_abellescampus <- data.frame(Shannon = shannondiv_abellescampus)
rownames(taula_shannon_abellescampus) <- rownames(data_shannon_abellescampus)
# Extreure dues noves columnes: tractament i data mostreig
taula_shannon_abellescampus <- taula_shannon_abellescampus %>%
mutate(Tractament = str_split(rownames(.), "_", simplify = TRUE)[, 2]) %>%
mutate(Data_Mostreig = str_extract(rownames(.), "(?<=-)[^-_]+-[0-9]+"))
print(taula_shannon_abellescampus)## Shannon Tractament Data_Mostreig
## Campus-juliol-2020_NS 2.611229 NS juliol-2020
## Campus-juliol-2020_S 1.905899 S juliol-2020
## Campus-Juliol-2021_NS 2.090031 NS Juliol-2021
## Campus-Juliol-2021_S 1.806531 S Juliol-2021
## Campus-Jun-2021_NS 2.701045 NS Jun-2021
## Campus-Jun-2021_S 2.080060 S Jun-2021
## Campus-maig-2020_NS 2.199104 NS maig-2020
# Fer el model
model_shannonabellescampus <- lmer(Shannon ~ Tractament + (1|Data_Mostreig), data = taula_shannon_abellescampus)
summary(model_shannonabellescampus)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Tractament + (1 | Data_Mostreig)
## Data: taula_shannon_abellescampus
##
## REML criterion at convergence: 1.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.03570 -0.47827 0.03651 0.57598 0.80378
##
## Random effects:
## Groups Name Variance Std.Dev.
## Data_Mostreig (Intercept) 0.03686 0.1920
## Residual 0.02461 0.1569
## Number of obs: 7, groups: Data_Mostreig, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.4004 0.1240 3.9717 19.363 4.42e-05 ***
## TractamentS -0.5097 0.1248 2.2656 -4.083 0.0444 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.398
TractamentS: quan el tractament és S, s’espera que l’índex de Shannon disminueixi 0,5097 respecte el tractament NS. És un resultat significatiu! (p-valor= 0.0444).
# Plot del model
plot(fitted(model_shannonabellescampus), resid(model_shannonabellescampus),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")MODELS LINEALS D’ANÀLISIS FUNCIONALS
Inter-tegular span (ITS) (model lineal mixt)
Abelles totals, amb Zona com a factor aleatori
#Obrir la base de dades amb la columna de les mesures de la ITS i filtrar dades per abelles
abelles_ITS <- dades_finals[dades_finals$Ordre == "HymenopteraAbella", ]
head(abelles_ITS)## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Model lineal mixt ITS (mida)
abelles_ITS_filtrades <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]
model_abelles_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|Zona), data = abelles_ITS_filtrades)
summary(model_abelles_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | Zona)
## Data: abelles_ITS_filtrades
##
## REML criterion at convergence: 1453.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3928 -0.4191 -0.2937 0.1489 5.4311
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 0.03073 0.1753
## Residual 0.48378 0.6955
## Number of obs: 683, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.78796 0.11569 2.44978 15.455 0.00167 **
## TractamentS -0.51521 0.05431 676.60529 -9.487 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.217
La ITS de les abelles disminueix -0.51521 quan el tractament és S (abelles més petites quan està segat). És significatiu (p-valor = 2e-16)
# Plot del model
# Prediccions del model sobre les dades:
predicted_values <- fitted(model_abelles_ITS)
observed_values <- abelles_ITS[!is.na(abelles_ITS$Inter_tegular_span),]$`Inter_tegular_span`
# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values, Predicted = predicted_values)
# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") + # Línea de referencia de 45 grados
labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
theme_minimal()Abelles del Campus, amb data_mostreig com a factor aleatori
# Filtrar les dades només per abelles del Campus
abellescampus_ITS <- abelles_ITS[abelles_ITS$Zona == "Campus",]
#Definir nova variable: "data_mostreig", combinació de mes i any de mostreig
abellescampus_ITS <- abellescampus_ITS %>%
mutate(data_mostreig = paste(Mes, Any, sep = "-"))
head(abellescampus_ITS)## # A tibble: 6 × 27
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 19 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>, data_mostreig <chr>
# Fer el model
model_abellescampus_ITS <- lmer(`Inter_tegular_span` ~ Tractament + (1|data_mostreig), data = abellescampus_ITS)
summary(model_abellescampus_ITS)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Inter_tegular_span ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## REML criterion at convergence: 1196.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3101 -0.5834 -0.2132 0.1382 5.4645
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.04629 0.2151
## Residual 0.42778 0.6540
## Number of obs: 594, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.54523 0.11430 3.29968 13.519 0.000529 ***
## TractamentS -0.36675 0.06534 537.59881 -5.613 3.19e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.243
Molt significatiu: la ITS de les abelles (mida) disminueix -0.36675 quan el tractament és S. Les abelles son més petites quan està segat.
# Fer el plot
# Prediccions del model sobre les dades:
predicted_values1 <- fitted(model_abellescampus_ITS)
observed_values1 <- abellescampus_ITS[!is.na(abellescampus_ITS$Inter_tegular_span),]$`Inter_tegular_span`
# Crear el data frame pel gràfic
plot_data <- data.frame(Observed = observed_values1, Predicted = predicted_values1)
# Graficar
ggplot(plot_data, aes(x = Predicted, y = Observed)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, color = "green", linetype = "dashed") +
labs(x = "Valors predits", y = "Valors observats", title = "Valors observats vs. predits") +
theme_minimal()ITS abelles per espècies totals:
# Crear nova base de dades: només una espècie per cada tractament.
especies_ITS <- abelles_ITS %>%
distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
head(especies_ITS)## # A tibble: 6 × 26
## Especimen Zona Mes `Mes num` Any Trampa Tractament Identificacio
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 6 Campus Juliol 7 2021 18 NS <NA>
## 2 8 Campus Juliol 7 2021 18 NS <NA>
## 3 9 Campus Juliol 7 2021 18 NS <NA>
## 4 11 Campus Juliol 7 2021 18 NS <NA>
## 5 27 Campus Juliol 7 2021 4 NS <NA>
## 6 159 Campus Juliol 7 2021 10 NS <NA>
## # ℹ 18 more variables: `Identificacio 2` <chr>, ID <chr>, Det. <chr>,
## # `M/F` <chr>, Genere <chr>, especie <chr>, Familia <chr>, Ordre <chr>,
## # NAME_JORDI <lgl>, Punxat <chr>, Observacions <chr>, ...20 <lgl>,
## # `nombre de etiqueta` <lgl>, ...22 <lgl>, Inter_tegular_span <dbl>,
## # Nesting_type <chr>, Lecty <chr>, Mostreig <chr>
# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especies_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especies_ITS)
summary(model_especies_ITS)##
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especies_ITS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3679 -0.7225 -0.3465 0.6686 3.5533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1199 0.1603 13.225 <2e-16 ***
## TractamentS -0.2894 0.2581 -1.121 0.266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.051 on 68 degrees of freedom
## (35 observations deleted due to missingness)
## Multiple R-squared: 0.01815, Adjusted R-squared: 0.003711
## F-statistic: 1.257 on 1 and 68 DF, p-value: 0.2662
No significatiu
ITS abelles per espècies del campus
# Crear nova base de dades: només una espècie per cada tractament.
especiescampus_ITS <- abellescampus_ITS %>%
distinct(ID, Tractament, .keep_all = TRUE)
# Dialectus -- dialectus
View(especiescampus_ITS)# Fer el model (lineal, no mixt: no podem considerar la zona com a factor aleatori ja que hem eliminat espècimens!)
model_especiescampus_ITS <- lm(`Inter_tegular_span` ~ Tractament, data = especiescampus_ITS)
summary(model_especiescampus_ITS)##
## Call:
## lm(formula = Inter_tegular_span ~ Tractament, data = especiescampus_ITS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2394 -0.6734 -0.3183 0.6081 2.8422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9914 0.1667 11.949 <2e-16 ***
## TractamentS -0.2301 0.2582 -0.891 0.377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9859 on 58 degrees of freedom
## (33 observations deleted due to missingness)
## Multiple R-squared: 0.01351, Adjusted R-squared: -0.003501
## F-statistic: 0.7941 on 1 and 58 DF, p-value: 0.3765
No significatiu.
Substrat de nidificació (GLMM binomial)
Nest_type en funció del tractament de les abelles totals (zona = factor aleatori)
abelles_ITS<-abelles_ITS[!is.na(abelles_ITS$Nesting_type),]
abelles_ITS["Nius"] = 0
abelles_ITS$Nius[abelles_ITS["Nesting_type"]=="S"]="S"
abelles_ITS$Nius[abelles_ITS["Nesting_type"]!="S"]="Not_S"
# Taula de contingència
taula_contingencia1 <- table(abelles_ITS$Tractament, abelles_ITS$Nius)
print(taula_contingencia1)##
## Not_S S
## NS 71 324
## S 14 291
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia1)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
ggtitle("Substrats de nidificació per tractament")# Model:
abelles_ITS$Nius <- factor(abelles_ITS$Nius)
# Nomes Zona o tambe mostreig, com a factors aleatoris
model_abelles_nest <- glmer(Nius ~ Tractament + (1|Zona) + (1|Mostreig), data = abelles_ITS, family = binomial)## boundary (singular) fit: see help('isSingular')
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Nius ~ Tractament + (1 | Zona) + (1 | Mostreig)
## Data: abelles_ITS
##
## AIC BIC logLik deviance df.resid
## 493.3 511.5 -242.6 485.3 696
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7820 0.2091 0.2282 0.4715 0.5146
##
## Random effects:
## Groups Name Variance Std.Dev.
## Mostreig (Intercept) 3.795e-02 0.1947978
## Zona (Intercept) 1.277e-10 0.0000113
## Number of obs: 700, groups: Mostreig, 6; Zona, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4843 0.1645 9.024 < 2e-16 ***
## TractamentS 1.6262 0.3450 4.714 2.43e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.442
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Significatiu! Les abelles del tractament S tenen més probabilitat (1.6262) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 2.43e-06).
# Plot del model
plot(fitted(model_abelles_nest), resid(model_abelles_nest),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Nest_type en funció del tractament de les abelles del Campus (data_mostreig = factor aleatori)
abellescampus_ITS<-abellescampus_ITS[!is.na(abellescampus_ITS$Nesting_type),]
abellescampus_ITS["Nius"] = 0
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]=="S"]="S"
abellescampus_ITS$Nius[abellescampus_ITS["Nesting_type"]!="S"]="Not_S"
# Taula de contingència
taula_contingencia2 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Nius)
print(taula_contingencia2)##
## Not_S S
## NS 63 291
## S 13 243
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia2)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Nesting_type") +
ggtitle("Substrats de nidificació per tractament")# Model:
abellescampus_ITS$Nesting_type <- factor(abellescampus_ITS$Nius)
model_abellescampus_nest <- glmer(Nesting_type ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)
summary(model_abellescampus_nest)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Nesting_type ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## AIC BIC logLik deviance df.resid
## 440.2 453.4 -217.1 434.2 607
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5162 0.2291 0.2360 0.4689 0.4997
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 0.02712 0.1647
## Number of obs: 610, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4968 0.1764 8.484 < 2e-16 ***
## TractamentS 1.5004 0.3860 3.887 0.000101 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.512
Significatiu! Les abelles del tractament S tenen més probabilitat (1.5004) de fer nius al sòl que de fer nius a altres llocs que no siguin el sòl. És un resultat significatiu! (p-valor = 0.000101).
# Plot del model
plot(fitted(model_abellescampus_nest), resid(model_abellescampus_nest),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Lèctia (GLMM binomial)
Lecty en funció del tractament de les abelles totals (Zona = factor aleatori)
# Taula de contingència
taula_contingencia3 <- table(abelles_ITS$Tractament, abelles_ITS$Lecty)
print(taula_contingencia3)##
## Oligolectic Polylectic
## NS 142 249
## S 38 263
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia3)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Lecty") +
ggtitle("Lecty per tractament")# Model:
abelles_ITS$Lecty <- factor(abelles_ITS$Lecty)
model_abelles_lecty <- glmer(Lecty ~ Tractament + (1|Zona), data = abelles_ITS, family = binomial)
summary(model_abelles_lecty)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Lecty ~ Tractament + (1 | Zona)
## Data: abelles_ITS
##
## AIC BIC logLik deviance df.resid
## 730.0 743.6 -362.0 724.0 689
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0678 -0.6616 0.3260 0.7371 1.5115
##
## Random effects:
## Groups Name Variance Std.Dev.
## Zona (Intercept) 0.7457 0.8635
## Number of obs: 692, groups: Zona, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2956 0.5350 0.553 0.581
## TractamentS 1.6319 0.2235 7.302 2.84e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.132
Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat molt significatiu (p-valor=), la qual cosa suggereix un fort efecte del tractament en la lecticitat de les abelles.
# Plot del model
plot(fitted(model_abelles_lecty), resid(model_abelles_lecty),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")Lecty en funció del tractament de les abelles del campus (data_mostreig = factor aleatori)
# Taula de contingència
taula_contingencia4 <- table(abellescampus_ITS$Tractament, abellescampus_ITS$Lecty)
print(taula_contingencia4)##
## Oligolectic Polylectic
## NS 127 223
## S 20 232
# Gràfic de barres:
df <- as.data.frame.table(taula_contingencia4)
ggplot(df, aes(x = Var1, y = Freq, fill = Var2)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Freq), position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Tractament", y = "Nº", fill = "Lecty") +
ggtitle("Lecty per tractament")# Model:
abellescampus_ITS$Lecty <- factor(abellescampus_ITS$Lecty)
model_abellescampus_lecty <- glmer(Lecty ~ Tractament + (1|data_mostreig), data = abellescampus_ITS, family = binomial)
summary(model_abellescampus_lecty)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Lecty ~ Tractament + (1 | data_mostreig)
## Data: abellescampus_ITS
##
## AIC BIC logLik deviance df.resid
## 468.8 482.0 -231.4 462.8 599
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2618 0.0518 0.1706 0.5796 1.2245
##
## Random effects:
## Groups Name Variance Std.Dev.
## data_mostreig (Intercept) 6.291 2.508
## Number of obs: 602, groups: data_mostreig, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2170 1.3182 1.682 0.0926 .
## TractamentS 0.6375 0.3271 1.949 0.0513 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TractamentS -0.091
Les abelles tenen una major probabilitat de ser polilèctiques quan estan sotmeses al tractament S. És un resultat marginalment significatiu (p-valor = 0.0513), la qual cosa suggereix un cert efecte del tractament en la lecticitat de les abelles.
# Plot del model
plot(fitted(model_abellescampus_lecty), resid(model_abellescampus_lecty),
xlab = "Valors ajustats", ylab = "Residus estandarditzats")